import re
import numpy as np
import matplotlib.pyplot as plt
import boto3
import rasterio
from rasterio.plot import show, show_hist
from skimage.exposure import equalize_adapthist
session = boto3.Session(profile_name='default')
scenes = [
'LC08_L1TP_205021_20180625_20180704_01_T1',
'LC08_L1TP_205021_20180727_20180731_01_T1'
]
def get_aws_band(scene, band, profile=False):
path, row = re.match(r'LC08_L1TP_(\d{3})(\d{3})', scene).groups()
prefix = f'c1/L8/{path}/{row}/{scene}'
with rasterio.Env(session=session):
with rasterio.open(
f's3://landsat-pds/{prefix}/{scene}_B{band}.TIF') as src:
if profile:
return src.profile
return src.read(out_shape=(src.height//10, src.width//10))
stacked_scenes = []
for scene in scenes:
red = get_aws_band(scene, 4)[0]
green = get_aws_band(scene, 3)[0]
blue = get_aws_band(scene, 2)[0]
stacked_scenes.append(np.stack([red, green, blue]))
We will use rasterio's plotting function (which wraps matplotlib)
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0], title="Scene 1", ax=axes[0])
show(stacked_scenes[1], title="Scene 2", ax=axes[1])
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers). Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55a42aeb8>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,65535)
axes[1].set_xlim(0,65535)
axes[0].set_ylim(0,140000)
axes[1].set_ylim(0,140000)
show_hist(
stacked_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Driver: GTiff/GeoTIFF
Files: data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Size is 8061, 8151
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-3],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32630"]]
Origin = (341085.000000000000000,6318015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 341085.000, 6318015.000) ( 5d36'53.37"W, 56d58'41.99"N)
Lower Left ( 341085.000, 6073485.000) ( 5d28'16.42"W, 54d47' 0.10"N)
Upper Right ( 582915.000, 6318015.000) ( 1d38' 6.39"W, 56d59'53.56"N)
Lower Right ( 582915.000, 6073485.000) ( 1d42'36.48"W, 54d48' 6.03"N)
Center ( 462000.000, 6195750.000) ( 3d36'28.13"W, 55d54'20.48"N)
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif | grep "Band"
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Most plotting tools work with 8 bit images (we have 16 bit):

fig, axes = plt.subplots(1,2, figsize=(30,15))
show(stacked_scenes[0].astype(np.uint8), title="Scene 1", ax=axes[0])
show(stacked_scenes[1].astype(np.uint8), title="Scene 2", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa569ab7128>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,50000)
axes[1].set_ylim(0,50000)
show_hist(
stacked_scenes[0].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
stacked_scenes[1].astype(np.uint8), bins=50, lw=0.0, stacked=False,
alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[1])
We have to be a bit smarter than that, we need to rescale!
def linear_rescale(image, in_range=[0, 65535], out_range=[1, 255]):
"""
Linear rescaling
"""
image = np.clip(image, in_range[0], in_range[1]) - in_range[0]
image = image / float(in_range[1] - in_range[0])
return (
image * (out_range[1] - out_range[0]) + out_range[0]
).astype(np.uint8)
linear_rescaled_images = [linear_rescale(stack) for stack in stacked_scenes]
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5539a90b8>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
OK, Landsat 8 doesn't actually make use of the full 16 bits, so we will rescale within the range we have reasonable values:
linear_rescaled_images = [
linear_rescale(stack, in_range=[0, 14000]) for stack in stacked_scenes
]
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5537ff438>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])

We can make better use of the colour range available, so lets stretch it.
def percentile_stretch(image, lower_perc=2, upper_perc=98):
mask = np.all(image != 0, axis=0).astype(np.uint8) * 255
image_out = np.zeros_like(image).astype(np.uint8)
for band in range(image.shape[0]):
input_range = np.percentile(
image[band][mask > 0],
(lower_perc, upper_perc)
).tolist()
image_out[band] = np.where(
mask,
linear_rescale(image[band],
in_range=input_range,
out_range=[0, 255]),
0
)
return image_out
percentile_stretched_images = [
percentile_stretch(stack) for stack in stacked_scenes
]
fig, axes = plt.subplots(1, 2, figsize=(30,15))
show(percentile_stretched_images[0], title="Scene 1", ax=axes[0])
show(percentile_stretched_images[1], title="Scene 2", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55394f828>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
percentile_stretched_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
percentile_stretched_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
That worked nicely with our clearish scene, but in the cloudy one, well the clouds look really interesting.
Contrast is lacking elsewhere...
Introducing Contrast Limited Adaptive histogram Equalization
![]()
def clahe_equalize(image, kernel_size=2000, clip_limit=0.03):
image_out = np.zeros_like(image).astype(np.float64)
for band in range(image.shape[0]):
image_out[band] = equalize_adapthist(
image[band], kernel_size=kernel_size, clip_limit=clip_limit
)
return image_out
raster_equalized_images = [
clahe_equalize(stack, clip_limit=0.05) for stack in stacked_scenes
]
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(
(raster_equalized_images[0]* 255).astype(np.uint8),
title="Scene 1", ax=axes[0]
)
show(
(raster_equalized_images[1]* 255).astype(np.uint8),
title="Scene 2", ax=axes[1]
)
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5536db048>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
(raster_equalized_images[0]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[0]
)
show_hist(
(raster_equalized_images[1]* 255).astype(np.uint8),
bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
title="Histogram", ax=axes[1]
)
The scenes still look very different, so not great for applying to all scenes, let's think about this in a physical sense.
We could try and remove the atmosphere, and ideally with an image based method (simpler, quicker, and no need for additional external data)
def find_dark_object_value(arr):
"""Find the value of the dark object
ie the first dark value that is not an outlier
"""
preval = None
step = arr.max() / 255.0
for val in np.unique(arr)[:100]:
if val == 0:
continue
if preval is not None and (val - preval) < step:
break
else:
preval = val
return preval
def dos(image):
image_out = np.zeros_like(image).astype(np.uint8)
# get the dark object for each band
for i, band in enumerate(image):
dark = find_dark_object_value(band)
image_out[i] = band - dark.astype(np.float)
image_out[image_out < 0] = 0
return image_out.astype(np.uint8)
prepped_scenes = [
linear_rescale(stack, [0, 14000]) for stack in stacked_scenes
]
dos_scenes = [dos(prepped_scene) for prepped_scene in prepped_scenes]
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55353c940>
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(linear_rescaled_images[1], title="Scene 2", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa568e31dd8>
fig, axes = plt.subplots(1,2, figsize=(30,15))
show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
<matplotlib.axes._subplots.AxesSubplot at 0x7fa568d85160>
fig, axes = plt.subplots(1,2, figsize=(30,10))
axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
dos_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
dos_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
histtype='stepfilled', title="Histogram", ax=axes[1])
def write_rendered_rgb(rgb, profile, outname_prefix):
"""
Write a rendered RGB to GeoTIFF.
"""
profile.update(dtype=rasterio.uint8, count=3, photometric='RGB')
outname = f'data/{outname_prefix}_rendered.tif'
with rasterio.open(outname, 'w', **profile) as dst:
for k, arr in enumerate(rgb):
dst.write_band(k+1, arr)
profile = get_aws_band(scenes[1], 4, True)
equalized_image_8bit = (raster_equalized_images[1] * 255).astype(np.uint8)
write_rendered_rgb(equalized_image_8bit, profile, "scene2")
!gdalinfo data/scene2_rendered.tif | grep Band
Band 1 Block=512x512 Type=Byte, ColorInterp=Red Band 2 Block=512x512 Type=Byte, ColorInterp=Green Band 3 Block=512x512 Type=Byte, ColorInterp=Blue